Théorie des Couches Minces Optiques

Ce document détaille la physique et les équations permettant de calculer les propriétés optiques d'un empilement de couches minces. Conçu pour les étudiants en physique numérique et optique, il tisse le lien entre l'intuition physique (interférences), les équations de Maxwell, et l'implémentation numérique via la méthode des matrices de transfert (TMM).

1. Introduction aux interférences en couches minces Niveau 1

Au lieu de sommer à la main la série infinie des réflexions multiples dans chaque couche (formule d'Airy), on utilise une approche globale dérivée de l'électromagnétisme : la méthode matricielle. Chaque couche est modélisée par une matrice 2×2 ; l'effet d'un empilement de 50 couches n'est que le produit de 50 matrices.

2. Conventions et notations Niveau 2

Pour garantir une lecture sans ambiguïté, l'ensemble des paramètres, variables et symboles utilisés dans ce document est regroupé ci-dessous selon leur domaine physique :

2.1 Paramètres de l'onde et grandeurs spatio-temporelles

2.2 Milieux, matériaux et empilement de couches

2.3 Électromagnétisme, champs et formalisme matriciel (Macleod)

2.4 Coefficients optiques, bilans énergétiques & photométrie

Convention de repère et de numérotation : l'axe $z$ est orienté vers le haut ($\uparrow$), dans le sens de croissance physique du dépôt. Le plan $z=0$ correspond à la surface du substrat, la couche 1 est collée au substrat ($0 \le z \le d_1$), et la couche $N$ est la dernière couche déposée, en contact avec l'air à $z = z_N = \sum_{j=1}^N d_j$.

3. Des équations de Maxwell à la matrice de transfert Niveau 2

Cette section établit rigoureusement, à partir des premiers principes de l'électromagnétisme, comment les équations de Maxwell se réduisent à l'équation de Helmholtz, puis comment les conditions aux limites et la condition d'onde sortante permettent de formuler la matrice de transfert $2 \times 2$ d'une couche mince.

3.1 Équations de Maxwell en régime harmonique

On adopte la convention temporelle de Macleod où tous les champs oscillent à la pulsation $\omega$ selon $e^{+i\omega t}$. Pour un champ électrique réel $\vec{\mathcal{E}}_{phys}(\vec{r}, t) = \mathrm{Re}[\vec{E}(\vec{r}) \, e^{i\omega t}]$, la dérivation temporelle devient un simple produit algébrique : $\frac{\partial}{\partial t} \longrightarrow +i\omega$.

Dans un milieu diélectrique ou conducteur linéaire, homogène, isotrope (LHI) et non magnétique ($\mu_r = 1$), sans charges libres macroscopiques ($\rho = 0$), les équations de Maxwell s'écrivent sous forme complexe :

La permittivité diélectrique complexe relative $\hat{\varepsilon}_r$ unifie la constante diélectrique $\varepsilon_r$ et la conductivité ohmique $\sigma$. Elle est directement reliée à l'indice de réfraction complexe $\hat{n} = n - i\kappa$ :

3.2 Équation d'onde et équation de Helmholtz

En appliquant l'opérateur rotationnel à l'équation de Maxwell-Faraday et en utilisant l'identité vectorielle $\nabla \times (\nabla \times \vec{E}) = \nabla(\nabla \cdot \vec{E}) - \Delta \vec{E}$ avec $\nabla \cdot \vec{E} = 0$ (milieu homogène), on obtient :

Comme la vitesse de la lumière dans le vide est $c = 1/\sqrt{\varepsilon_0 \mu_0}$ et le nombre d'onde dans le vide est $k_0 = \omega/c = 2\pi/\lambda$, on aboutit à l'équation vectorielle de Helmholtz :

Considérons une onde plane incidente se propageant dans le plan d'incidence $(x,z)$, où $z$ est l'axe perpendiculaire aux couches et $x$ l'axe parallèle aux interfaces. Par invariance par translation selon l'axe transverse $y$, on a $\partial/\partial y = 0$. La dépendance spatiale selon $x$ est commune à toutes les couches et fixée par le milieu incident :

L'équation de Helmholtz se réduit alors à une équation différentielle scalaire unidimensionnelle selon la coordonnée $z$ pour chaque composante :

3.3 Conditions de continuité aux interfaces (Conditions aux limites)

À chaque interface $z = z_k$ séparant deux milieux d'indices $\hat{n}_1$ et $\hat{n}_2$, les équations de Maxwell imposent des relations de continuité strictes. En intégrant Maxwell-Faraday et Maxwell-Ampère sur un contour rectangulaire infinitésimal chevauchant l'interface (largeur $\Delta l$, hauteur $\Delta h \to 0$), le flux magnétique et le courant traversant la boucle tendent vers zéro :

En l'absence de charges et de courants superficiels libres, les composantes tangentielles du champ électrique et du champ magnétique sont strictement continues :

3.4 Ondes progressives et Condition d'onde sortante (Sommerfeld)

Dans chaque couche homogène $j$ (délimitée entre sa base $z_{j-1}$ et son sommet $z_j = z_{j-1} + d_j$), la solution générale de l'équation de Helmholtz 1D est la superposition exacte d'une onde se propageant vers le bas (vers le substrat, dans le sens des $z$ décroissants, notée $E_j^-$) et d'une onde réfléchie se propageant vers le haut (vers l'air, dans le sens des $z$ croissants, notée $E_j^+$) :

3.5 Dérivation analytique de la matrice caractéristique 2×2

Pour relier le champ magnétique tangentiel $H_\parallel$ au champ électrique tangentiel $E_\parallel$, on utilise l'équation de Maxwell-Faraday $\nabla \times \vec{E} = -i\omega\mu_0 \vec{H}$ :

• Polarisation TE (s) : $\vec{E} = (0, E_y, 0)$. La composante $x$ de Maxwell-Faraday donne $-\partial_z E_y = -i\omega\mu_0 H_x$, donc $H_x = \frac{1}{i\omega\mu_0} \frac{\partial E_y}{\partial z} = \eta_s (E^- e^{+ik_z z} - E^+ e^{-ik_z z})$ avec l'admittance $\eta_s = \hat{n}\cos\theta$.
• Polarisation TM (p) : $\vec{H} = (0, H_y, 0)$. La composante $x$ de Maxwell-Ampère donne $-\partial_z H_y = i\omega\varepsilon_0\hat{n}^2 E_x$, donc $E_x = -\frac{1}{i\omega\varepsilon_0\hat{n}^2} \frac{\partial H_y}{\partial z}$ avec l'admittance $\eta_p = \frac{\hat{n}}{\cos\theta}$.
• Attention aux manuels : ces signes sont ceux de la convention $e^{+i\omega t}$ du § 3.1. Avec $e^{-i\omega t}$ (Born & Wolf, Jackson), les deux signes s'inversent ; le résultat final, lui, est le même.

Considérons une couche $j$ d'épaisseur $d$ comprise entre son interface supérieure d'entrée $a$ ($z = z_j$) et son interface inférieure de sortie $b$ ($z = z_{j-1} = z_j - d$) :

1. À l'interface supérieure d'entrée ($z = z_j$), en rapportant cette fois les deux amplitudes $E^\pm$ à cette interface (au § 3.4, $E_j^+$ était rapportée à la base de la couche : simple choix d'origine des phases, sans effet sur la matrice) :

Ce qui permet d'exprimer les amplitudes d'ondes en fonction des champs totaux d'entrée :

2. À l'interface inférieure de sortie ($z = z_{j-1} = z_j - d$), en posant le déphasage $\delta = k_z d$ :

3. Inversion matricielle pour exprimer les champs d'entrée $(E_a, H_a)$ en fonction des champs de sortie $(E_b, H_b)$ :

4. Méthode des matrices de transfert (TMM) Niveau 2

4.1 Matrice caractéristique d'une couche $j$

Pour une couche d'indice complexe $\hat{n}_j$ et d'épaisseur $d_j$, la matrice $M_j$ s'écrit :

Où $\delta_j$ est la phase optique accumulée : $\delta_j = \frac{2\pi}{\lambda} (\hat{n}_j) d_j \cos \theta_j$.

if n_cos_th.imag > 0: 
    n_cos_th = -n_cos_th  # Force l'atténuation physique

4.2 Matrice globale de l'empilement

La matrice globale $M$ du système est le produit matriciel des matrices de chaque couche. On multiplie en partant de la base sur le substrat ($z=0$, couche 1) vers le sommet en contact avec l'air ($z=z_N$, couche $N$) :

5. Polarisation et admittance optique Niveau 2

L'admittance optique $\eta_j$ caractérise le rapport entre les champs tangentiels ($H_\parallel / E_\parallel$). En incidence oblique, elle dépend de la polarisation de la lumière :

6. Coefficients de Fresnel et puissances (R, T) Niveau 2

Pour l'empilement complet, on définit le vecteur caractéristique $\begin{pmatrix} B \\ C \end{pmatrix}$ en multipliant la matrice globale $M$ par le vecteur des conditions aux limites du substrat :

Les coefficients d'amplitude en réflexion ($r$) et en transmission ($t$) s'écrivent :

On en déduit les grandeurs mesurables expérimentalement, à savoir la réflectance $R$ et la transmittance $T$ (avec conservation de l'énergie $R+T+A=1$, où $A$ est l'absorption) :

7. Substrat d'épaisseur finie Niveau 2

Les équations précédentes supposent un substrat de sortie semi-infini. En réalité, une lame de verre possède une face arrière. En raison de l'épaisseur importante du substrat (millimétrique, $\gg \lambda$), la cohérence optique spatiale est perdue. On ne somme donc plus les amplitudes complexes (interférences), mais les intensités des rebonds multiples incohérents dans le verre.

À retenir : en optique linéaire passive, la transmittance de puissance est réciproque. Autrement dit, si elle est bien définie et bien normalisée, on a $T'_{front} = T_{front}$. Ce qui peut changer avec le sens d'illumination, ce sont surtout $R$ et $A$.

Important : dans le cas d'un substrat fini, $T_{total}$ n'est pas la même chose que $T_{front}$. Il faut appliquer la formule complète ci-dessus ; sinon, on commet une erreur de physique.

Exemple rapide : un verre nu réfléchit environ 4% par face ($R_{front} = R'_{front} = R_{back} \approx 0.04$). La réflectance totale d'une vitre simple vaut alors environ $7.7\%$. La transmittance totale doit, elle aussi, être recalculée avec la formule de substrat fini.

8. [Avancé] Profil d'intensité du champ électrique $|E(z)|^2$ Niveau 3

Calculer le profil d'intensité à l'intérieur de l'empilement est fondamental, notamment pour évaluer les seuils de dommage laser (LIDT) ou optimiser l'absorption dans des cellules solaires.

Méthode : On connaît le champ transmis dans le substrat à $z=0$ (qui vaut l'amplitude transmise $t$). On peut remonter mathématiquement le long de l'axe $z$ (dans le sens de croissance des couches, de $z=0$ jusqu'au sommet $z=z_N$) couche après couche. Dans la couche $j$ ($z_{j-1} \le z \le z_j$), on part des champs à sa base, obtenus en appliquant successivement $M_1, \dots, M_{j-1}$ au vecteur de départ $t\,(1, \eta_{sub})^T$ :

Où $\delta_j(z) = \frac{2\pi}{\lambda} \hat{n}_j (z - z_{j-1}) \cos \theta_j$ est le déphasage accumulé depuis la base de la couche $j$. Pour une couche unique ($j = 1$, $z_0 = 0$), le départ est directement le champ au substrat, $E(0) = t$.

9. Exemples d'application numérique Niveau 1 & 2

Cas 1 : Couche Anti-Reflet (AR) quart-d'onde ($\lambda/4$)

L'objectif est d'annuler les $\sim 4\%$ de réflexion parasite du verre. On dépose une épaisseur telle que le déphasage aller-retour soit de $\pi$ (interférence totalement destructive). L'épaisseur optique requise est $nd = \lambda / 4$.

La réflexion s'annule complètement si l'indice de la couche vérifie $n_{couche} = \sqrt{n_{air} \cdot n_{verre}}$. Pour un verre standard, on voudrait $n \approx 1.22$. Le Fluorure de Magnésium (MgF₂, $n=1.38$) est le matériau approchant le plus couramment utilisé en optique classique pour s'approcher de cet idéal.

Cas 2 : Filtre de validation à 5 couches (avec métal absorbant)

Ce cas de test complexe permet de valider le code Python, car il inclut des angles obliques, de la polarisation et de l'absorption. Il sert surtout à vérifier deux choses :

• la conservation énergétique $R + T + A = 1$,
• la réciprocité de la transmittance de puissance : $T = T'$, même dans un milieu absorbant, dès lors que la comparaison est faite avec la même normalisation physique.

• Incidence : Air ($n=1.0$) | Substrat : Verre ($n=1.52$) à $\lambda$ = 550 nm.
• Empilement (numéroté du substrat vers l'air) :
  1. TiO₂ ($n=2.35$), 85 nm
  2. SiO₂ ($n=1.46$), 110 nm
  3. Ag ($n=0.12 - 3.45i$), 15 nm (Couche de métal absorbant)
  4. TiO₂ ($n=2.35$), 85 nm
  5. SiO₂ ($n=1.46$), 90 nm

10. Implémentation Numérique (Python) Niveau 2

Le code ci-dessous résume mathématiquement la théorie exposée. Il prend en compte l'absorption, l'angle d'incidence, les polarisations et corrige le problème de la racine complexe pour les ondes évanescentes (voir § 4.1).

import numpy as np

def compute_tmm(n_inc, n_exit, n_layers, d_layers, wl, theta0=0.0, pol='s'):
    """
    Calcule les propriétés optiques d'un empilement de couches minces (TMM).
    
    Arguments :
      n_inc, n_exit : Indices complexes du milieu incident (air) et du substrat.
      n_layers      : Liste des indices complexes des couches (n - i*k), de la
                      couche 1 (contre le substrat) à la couche N (contre l'air).
      d_layers      : Liste des épaisseurs physiques des couches (en nm), même ordre.
      wl            : Longueur d'onde de travail (en nm).
      theta0        : Angle d'incidence dans le milieu incident (en radians).
      pol           : Polarisation ('s' pour TE, 'p' pour TM).
      
    Retourne :
      R (Réflectance), T (Transmittance), A (Absorption)
    """
    
    # =====================================================================
    # Étape 1 : Protection de la Convention de Macleod (Assure n - ik)
    # =====================================================================
    # En optique des couches minces, l'onde s'écrit formellement exp(i(wt - kz)).
    # Pour que le milieu absorbe l'énergie (et ne l'amplifie pas), la partie 
    # imaginaire de l'indice complexe doit être strictement négative.
    def enforce_macleod(n): 
        return n if n.imag <= 0 else n.real - 1j * n.imag
    
    n_inc = enforce_macleod(n_inc)
    n_exit = enforce_macleod(n_exit)
    n_layers = [enforce_macleod(n) for n in n_layers]
    
    # Nombre d'onde dans le vide
    k0 = 2 * np.pi / wl
    
    # =====================================================================
    # Étape 2 : Loi de Snell-Descartes et Angles Complexes (Correction TIR)
    # =====================================================================
    # L'invariant de Snell-Descartes (n * sin(theta)) est conservé partout.
    S0 = n_inc * np.sin(theta0) 
    
    def safe_cos_th(n_val):
        # On utilise l'identité trigonométrique classique : 
        # cos(theta) = sqrt(n^2 - (n*sin(theta))^2) / n
        n_cos_th = np.sqrt(n_val**2 - S0**2 + 0j)
        
        # ATTENTION AU PIÈGE NUMÉRIQUE (Coupure complexe / Branch cut) :
        # En cas de Réflexion Totale Interne (TIR) ou d'onde évanescente, 
        # numpy peut renvoyer une partie imaginaire positive, ce qui simule 
        # un gain non-physique d'énergie. On force la décroissance physique :
        if n_cos_th.imag > 0: 
            n_cos_th = -n_cos_th 
        return n_cos_th / n_val

    # Calcul des cosinus (complexes si TIR/absorption) pour chaque milieu
    cos_inc = np.cos(theta0) 
    cos_exit = safe_cos_th(n_exit)
    cos_layers = [safe_cos_th(n) for n in n_layers]
    
    # =====================================================================
    # Étape 3 : Calcul des Admittances optiques (eta)
    # =====================================================================
    # L'admittance relie le champ magnétique au champ électrique (H = eta * E)
    def get_eta(n_val, cos_th, pol):
        if pol.lower() == 's':
            # Polarisation TE (Champ E transverse) : projection par multiplication
            return n_val * cos_th       
        else:
            # Polarisation TM (Champ H transverse) : projection par division
            return n_val / cos_th       
            
    eta_inc = get_eta(n_inc, cos_inc, pol)
    eta_exit = get_eta(n_exit, cos_exit, pol)
    eta_layers = [get_eta(n_layers[i], cos_layers[i], pol) for i in range(len(n_layers))]
    
    # =====================================================================
    # Étape 4 : Multiplication Matricielle (Du substrat vers l'air)
    # =====================================================================
    # On initialise la matrice système globale à la matrice Identité
    M_sys = np.eye(2, dtype=complex)
        
    for i in range(len(n_layers)):
        # Calcul du déphasage optique local (delta) accumulé dans la couche i
        delta = k0 * n_layers[i] * d_layers[i] * cos_layers[i]
        c, s = np.cos(delta), 1j * np.sin(delta)
        eta = eta_layers[i]
        
        # Construction de la matrice caractéristique 2x2 de la couche locale i
        M_i = np.array([
            [c, s / eta], 
            [s * eta, c]
        ], dtype=complex)
        
        # Accumulation : M_sys = M_i * M_sys
        # L'ordre de multiplication est crucial : on part du substrat (sortie) 
        # vers le milieu incident (entrée), donc de droite à gauche.
        M_sys = np.dot(M_i, M_sys)

    # =====================================================================
    # Étape 5 : Extraction des coefficients d'amplitude de Fresnel
    # =====================================================================
    # On applique la matrice globale aux conditions aux limites du substrat
    # Équation matricielle : [B, C]^T = M_sys * [1, eta_exit]^T
    B = M_sys[0,0]*1 + M_sys[0,1]*eta_exit
    C = M_sys[1,0]*1 + M_sys[1,1]*eta_exit
    
    # Calcul des coefficients complexes de réflexion (r) et transmission (t)
    r = (eta_inc * B - C) / (eta_inc * B + C)
    t = (2.0 * eta_inc) / (eta_inc * B + C)
    
    # =====================================================================
    # Étape 6 : Bilan des puissances (Conservation de l'énergie)
    # =====================================================================
    # La réflectance se mesure dans le même milieu, donc un simple module carré
    R = abs(r)**2
    
    # Correction du vecteur de Poynting (Flux d'énergie) :
    # La transmittance nécessite de comparer la partie réelle des admittances
    # des milieux extrêmes car l'onde change de milieu de propagation.
    T_factor = (eta_exit.real) / (eta_inc.real) if eta_inc.real > 1e-12 else 0.0
    T = T_factor * abs(t)**2
    
    # L'énergie manquante est nécessairement absorbée par l'empilement
    A = 1.0 - R - T
    
    return R, T, A

# --- Exemple d'utilisation simple (Verre nu) ---
# Pour un verre d'indice 1.5 dans l'air (1.0) à 550 nm en incidence normale :
# R, T, A = compute_tmm(1.0, 1.5, [], [], 550) 
# print(f"Réflectance du verre nu : {R*100:.2f}%") # Affichera ~4.00%